The question is how much G x E we see in the RNAseq data.
To get started on this I ran the file “ProcessCount.R”. For this the data was TMM normalized, voom transformed, and then a gtXenvironment model was fit in limma.
library(limma)
library(tidyverse)
[30m── [1mAttaching packages[22m ─────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✔[30m [34mggplot2[30m 2.2.1 [32m✔[30m [34mpurrr [30m 0.2.4
[32m✔[30m [34mtibble [30m 1.3.4 [32m✔[30m [34mdplyr [30m 0.7.4
[32m✔[30m [34mtidyr [30m 0.7.2 [32m✔[30m [34mstringr[30m 1.2.0
[32m✔[30m [34mreadr [30m 1.1.1 [32m✔[30m [34mforcats[30m 0.2.0[39m
[30m── [1mConflicts[22m ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31m✖[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(stringr)
library(modelr)
library(broom)
Attaching package: ‘broom’
The following object is masked from ‘package:modelr’:
bootstrap
library(multidplyr)
Get the limma objects
load("voom.int.Rdata")
How many genes show GxE?
generate vector of interaction coefficients and then topTable
colnames(fit.int$coefficients)
int.coef.n <- colnames(fit.int$coefficients) %>% str_which(":")
topTable(fit.int,coef = int.coef.n)
Not working
let’s just do an ANOVA for each gene
counts <- as.tibble(counts.voom.gr$E) %>%
mutate(geneIndex=1:n())
counts[1:10,1:10]
counts.long <- counts %>% gather(key="name",value="expression", -geneIndex) %>%
mutate(gt=factor(str_extract(name,"RIL_[0-9]*")),
env=factor(str_extract(name,"CR|UN")),
rep = factor(str_extract(name,"Rep[0-9]"))) %>%
select(-name) %>%
group_by(geneIndex) %>% nest()
counts.long